Impact of the Grain for Green Project on water resources and ecological water stress in the Yanhe River Basin

The Grain for Green project (GGP), initialized by the Chinese government in 1999, has achieved substantial achievements accompanied by a decrease in surface runoff on the Loess Plateau, but the impacts of large-scale afforestation on regional water resources are uncertain. Hence, the objective of this study was to explore the impact of land use change on generalized water resources and ecological water stress using the blue and green water concepts, taking the Yanhe River Basin as the case study. The Soil and Water Assessment Tool (SWAT) was applied to quantify the green water and blue water, which are defined as generalized water resources. The ecological water requirement of vegetation (forest and grass), agricultural water footprint and virtual water flow are considered regional water requirements. The land use types of 1980 (Scenario I) and 2017 (Scenario II) were entered into the SWAT model while keeping the other parameters constant to isolate the influence of land use changes. The results show that the average annual differences in blue, green and generalized water resources were -72.08 million m3, 24.34 million m3, and -47.74 million m3, respectively, when the simulation results of Scenario II were subtracted from those of Scenario I, which shows that land use change caused by the GGP led to a decrease in blue and generalized water resources and an increase in green water resources. Surface runoff in Scenario I was more than that in Scenario II in all of the years of the study period from 1980–2017, and green water storage in Scenario I was more than that in Scenario II in all of the years of the study period except in 1998; although lateral flow in Scenario I was less than that in Scenario II except in 2000 and 2015, as was groundwater runoff in 1992, 2000 and 2015, and green water flow in 1998. Blue water flow, green water storage and green water flow in Scenario II were less than those in Scenario I in the whole basin, 12.89 percent of the basin and 99.21 percent of the basin, respectively. The total water footprint increased from 1995 to 2010 because the forest water footprint increased significantly in this period, although the agricultural water footprint and grass water footprint decreased. The ecological water stress index values had no obvious temporal change trends in either land use scenario, but the ecological water stress index in Scenario II was greater than that in Scenario I, which illustrates that the GGP led to an increase in ecological water stress from the perspective of generalized water resources.


Introduction
The Chinese government launched the Grain for Green project (GGP) in the 1990s to control soil erosion and water loss in the Loess Plateau [1][2][3]. After more than two decades of vegetation restoration, soil erosion caused by unreasonable land use has been curbed, and the ecological environment of the region have been greatly improved [4,5]. Theoretically, vegetation restoration can enhance vegetation coverage, increase precipitation interception and water retention, decrease soil erosion and thus improve ecosystem stability. However, forests consume more water than other vegetation types, such as agricultural crops and natural grasslands [6]; therefore, accompanied by the enhancement of vegetation coverage since 2000 in China, runoff has shown a significant decrease in Haihe, Liaohe, Songhua Jiang, Hanjiang and the Yellow River [7][8][9]. In particular, runoff in the middle reaches of the Yellow River, which has the most obvious vegetation restoration achievements, has decreased sharply, and the annual flow at the Huayuankou station decreased from 55.9 billion m 3 (1970s) to 45.2 billion m 3 (2010)(2011)(2012)(2013)(2014)(2015) [10]. Some results indicate that large-scale vegetation restoration on the Loess Plateau has positive impacts on soil erosion control and the ecological environment and negative impacts on streamflow [11][12][13]. Studies in other areas of the world have also demonstrated that an increase in vegetation coverage will lead to higher interception loss, which is the main reason for streamflow reductions [14,15]. Some studies have also indicated that unsuitable vegetation types and overlooking biodiversity will bring about soil desiccation on the Loess Plateau [16,17]. Studies have also presented the negative effects of afforestation on underground water resources [18] and ecological water deficits because of afforestation [19]. Why do streamflow decrease and where does the water go? Has the soil become dryer and has water stress become more serious because of the GGP? These are important problems need to be discussed.
The blue and green water concepts proposed by Falkenmark [20] provide new theories and ideas for water resource management, especially in arid and semiarid regions. A large amount of blue water converting to green water is one of the important causes of the streamflow decline on the Loess Plateau [21]. Vegetation coverage improvement significantly reduces streamflow, and vegetation restoration is close to the threshold of the water resources carrying capacity from the perspective of blue water [22]. However, it can reduce the water in sediment transport, and the virtual water embodied in green plants is far greater than streamflow reduction from the perspective of green water [23]. The water footprint (WF) proposed by Hoekstra and Hung [24] represents direct and indirect measurements of water appropriation by human beings. It quantifies blue and green water consumption in a river basin or a specified region and is a new approach to assess the sustainable water use of economic production sectors or regions [25][26][27].
The Yanhe River Basin in the Loess Plateau is the first tributary of the Yellow River and is in a semi-arid area having serious water scarcity and severe soil erosion. This watershed was one of the earliest and fastest areas in the whole country to convert cultivated land (grasslands) to forests and the vegetation restoration effect has been significant since the implementation of the GGP. During the period of 2000-2017, forests increased by 2357.6 km 2 , cultivated land decreased by 2116.2 km 2 , urban land increased by 222.1 km 2 , and water, grassland and other grass, urban use land, waters and barren land. Digital evaluation model (DEM) data with 30 m resolution were provided by the Geospatial Data Cloud (http://www.gscloud.cn/). The daily meteorological data were extracted from China's meteorological data sharing service system (http:// cdc.cma.gov.cn/home.do). Soil data with an accuracy of 1:1000000 were obtained from the Data Center for Resources and Environmental Sciences (http://vdb3.soil.csdb.cn/). Annual and monthly streamflow data of the hydrological stations were provided by the Loess Plateau SubCenter, National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://loess.geodata.cn). Various crop yields, populations, and grain consumption were obtained from Yan'an Statistical Yearbooks , which were published by China Statistics Press; forest and grassland areas were extracted from land use maps. All data in the manuscript were collected by the authors from the public information website and data center, the website was attached to the data statement. We confirmed that others would be able to access these data in the same manner as the authors.
land cover types did not change remarkably. Numerous studies have attempted to evaluate the impact of vegetation restoration on water resources and have focused on improving the watershed ecological environment, ameliorating soil properties, and changing streamflow and sediment [28][29][30]. However, there has been almost no study exploring the impact of vegetation restoration on water stress from the aspect of the water footprint. The WF can assess natural water resource availability, support optimal allocation among different regions and improve watershed sustainability. It can evaluate water demand and consumption more comprehensively from the perspective of blue water, green water and WF. Therefore, the objectives of this study were to (1) analyze the spatial-temporal characteristics of land use change during the period of 1980-2017 in the Yanhe River Basin and describe GGP achievements; (2) quantify water balance elements and analyze the spatial-temporal characteristics of green and blue water over the whole basin and subbasins based on calibrated and validated SWAT model simulation results; (3) investigate the temporal characteristics of agricultural WF, ecological WF and virtual water flow; and (4) calculate the ecological water stress index based on generalized water resources and the water footprint and probe the impact of GGP on regional water stress.

Study area
The Yanhe River Basin (36.21'-37.19'N, 108.38'-110.29'E) is located on the Loess Plateau in northern Shaanxi Province in China and is a first-order tributary of the Yellow River (Fig 1). With an area of 7785 km 2 , the watershed has a warm temperate continental semiarid monsoon climate. The annual mean precipitation of the watershed is 520 mm, with 75% being concentrated from June to September, and the mean annual temperature varies from 8.8 to 10.2 degrees Celsius [31]. The yearly streamflow at Gan'guyi was 220 million m 3 , and the major land use and land cover types of the watershed were forestland, shrub land, grassland, cultivated land, construction land, water body, and bare land. The annual water resource per capita value in the basin was 375 m 3 , accounting for 28 percent of that of Shaanxi Province and 17 percent of that of China. This river has a large sediment content and serious point source and nonpoint source pollution. The water resources in the watershed are in acute shortage, and the ecological environment is fragile.

Data sources
The land uses of the region in 1980 (before the GGP) were obtained from the Resource and Environment Science and Data Center (http://www.resdc.cn), and the land uses of 2017 (after the GGP) were interpreted from Landsat OLI images. The image data having 30 m resolution and zero cloud cover were collected from the Geospatial Data Cloud (http://www.gscloud.cn/). The land uses of 1990, 2000, 2005, and 2010 were also from the Resource and Environment Science and Data Center (http://www.resdc.cn). The land use types in the Yanhe River basin are classified into cultivated land, forest, grass, urban use land, waters and barren land. Digital evaluation model (DEM) data with 30 m resolution were provided by the Geospatial Data Cloud (http://www.gscloud.cn/). The daily meteorological data were extracted from China's meteorological data sharing service system (http://cdc.cma.gov.cn/home.do). Soil data with an accuracy of 1:1000000 were obtained from the Data Center for Resources and Environmental Sciences (http://vdb3.soil.csdb.cn/). Annual and monthly streamflow data of the hydrological stations were provided by the Loess Plateau SubCenter, National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://loess.geodata.cn). Various crop yields, populations, and grain consumption were obtained from Yan'an Statistical Yearbooks (1995-2018), which were published by China Statistics Press; forest and grassland areas were extracted from land use maps.

SWAT model
SWAT is a physically based, semidistributed model that has the advantage of simulating the quality and quantity of both surface and ground water as well as predicting the impact of land cover change, land management practices and climate change [32][33][34]. The model has been successfully applied to calculate water yield and evaluate water quality in small watersheds and large river basins [35,36]. Simulation results can be applied to quantify the spatial and temporal characteristics of blue and green water in different parts of the world, such as the Savannah  The sequential uncertainty fitting (SUFI-2) algorithm (Abbaspour, 2015) [40,41] in SWAT calibration and uncertainty programs (SWAT-CUP) were used to perform parameter estimation and sensitivity analysis. The algorithm describes all parameters and tries to simulate the observed data within the 95% confidence interval (95 PPU) by using an iterative process. The primary target of calibration is to identify the sensitive parameters in the watershed that control surface runoff. The P-factor and R-factor are used to quantify the uncertainty in this algorithm. The P-factor measures the percentage of the measured data with ninety-five percent prediction certainty. The R-factor represents the average distance of 95 PPU divided by the standard deviation of the measured data. P-factor values range between 0 and 1, and the R-factor values range from 0 to infinity; the closer the P-factor is to 1 and the closer the R-factor is to 0, the better the calibration and uncertainty analysis results [42]. R 2 (the coefficient of determination) and NSE (the Nash-Sutcliffe efficiency) are used to evaluate the accuracy of the SWAT model. R 2 values range from 0 to 1, and if the value is close to 1, it indicates that the simulated data are close to the observed data and vice versa [43].

Estimating blue and green water resources
Blue water resources in each subbasin are equal to water yield (WYLD), which is the sum of surface runoff (SURQ), lateral flows (LATQ) and ground water recharge (GWQ) [44]. Green water resources consist of green water flow (actual evapotranspiration, ET) and green water storage (soil water content, SW). The annual SURQ, LATQ, GWQ, ET and SW were obtained from simulation results using the calibrated and validated SWAT model [45,46].
To investigate the influences of land use change caused by the GGP on blue, green water and each hydrological element, two scenarios are set up by changing land use types while keeping the other input parameters in SWAT model unchanged (scenario I: land uses of 1980; scenario II: land uses of 2017). Only the impact of land use changes was considered when the simulation results of scenario II were subtracted from those of scenario I.

WF
The WF concept quantifying the volumetric freshwater consumption of products is distinguished as green, blue and gray WFs [20]. Green WF (WFgreen) is defined as rainwater that is stored in the soil and evaporated or consumed during production. Blue WF (WFblue) refers to surface and groundwater that are consumed or evaporated during production. The gray water footprint (WFgray) expresses the volume of water needed to dilute pollutants to achieve allowed values in the receiving water bodies. The sum of WFgreen and WFblue demonstrates the quantity of total freshwater consumption during production, while WFgray indicates the degradation of water quality. The goal of this study has been to explore water quantity changes caused by vegetation restoration, so WFgray is not considered during the calculation process.
Agricultural crop WF. The crop WF is the sum of WF green and WF blue , which are determined by the following equations [47]: where WF green and WF blue are the green WF (m 3 /t) and blue WF (m 3 /t) during the crop growth season, respectively; CWU green and CWU blue are the green and blue water uses (m 3 /ha); 10 is a constant to convert the water depth (mm) to the water volume (m 3 /ha); Y is the crop yield (m 3 /ha); and ET green and ET blue are defined as the evaporative demand satisfied by green and blue water, respectively. ET green and ET blue are calculated as [48]: where ET c represents the actual evapotranspiration of crops from sowing day to the harvest and P eff is the effective precipitation [49].
where K c is the crop coefficient, ET 0 is the potential reference crop evapotranspiration and is calculated by the Penman-Monteith formula, and P is the precipitation of ten days. Agricultural virtual water flow calculation. Due to the lack of statistics on regional trade and storage of agricultural products in China, the calculation of agricultural virtual water flow in this study was based on the assumption that the regional agricultural product storage remains unchanged and the water footprints per unit mass of local agricultural product inputs and outputs are equal. The equation is as follows [50]: where W v,f is the regional virtual flow; while W v,pro , and W v,con are the amounts of virtual water production and water consumption, respectively. W v,pro and W v,con can be calculated by multiplying the agricultural virtual water per unit mass times the amounts of grain production and grain consumption, respectively. If W v,f is less than 0, it indicates that the region inputs virtual water from outside; if W v,f is greater than 0, it indicates that the region exports virtual water, which may aggravate regional water resources pressure; if is W v,f is equal to 0, it indicates that there is no transfer of agricultural products between the region and the outside. Vegetation WF. The equation of vegetation WF is [51]: where WF green is the vegetation WF (m 3 ); A p is the area of vegetation coverage (m 2 ); ET p is the vegetation evapotranspiration (mm/day) under restricted circumstances; and p is the vegetation type. Vegetation evapotranspiration is less than the potential evapotranspiration when the soil water content is below a specified threshold, and the effect is determined by the soil moisture limitation coefficient (K s ). Thus, the calculation equation of ET p is as follows [52]: where ET 0 is the potential reference evapotranspiration, K c is the vegetation water demand coefficient, K s is the soil moisture limitation coefficient, and j is the soil type.
Based on previous studies on the Loess Plateau [53,54], the K c values of forest and grassland are 0.765 and 0.65, respectively. The Yanhe River Basin contains silty soil and sandy loam soil, and the K s values of the two soil types are 0.537 and 0.556, respectively.

Ecological water stress index (EWSI)
Raskin [55] proposed a criterion using the ratio between water demand and available water resources to estimate water scarcity, and it has been widely used to evaluate global and regional water resources [56,57]. The water stress index (WSI) can provide information on the management of freshwater resources [58]. In this paper, the EWSI is calculated as the ratio of the ecological water footprint to generalized water resources (Fig 2). Generalized water resources are the sum of blue water and green water simulated by the SWAT model, and the ecological water footprint, which contains the agricultural WF, forest WF, grass WF and virtual water flow.

Land use change in Yanhe River Basin
The Yanhe River Basin land use types of 1980 and 2017 are shown in Table 1   The diagonal value in the land use conversion matrix (Table 2) is the unchanged area of each land use type. 5132.30 km 2 of land use type has changed during the period from 1980 to 2017, 2937.70 km 2 of cultivated land has changed into other land use types and the area proportion is 57.24%; 1198.30 km 2 has been converted into forest and 1586.44 km 2 has been converted into grassland. Grassland is the second land use conversion type (1843.56 km 2 ) and it was mainly converted into forest (1337.85 km 2 ), cultivated land (388.01 km 2 ) and urban use land (103.85 km 2 ). The transfer-out land use type was cultivated land and grassland, and the transfer proportion accounted for as high as 93.16% of the whole basin.
A total of 2545.06 km 2 of forest was developed from other land use types since 1980, and most of the conversion was from grass land and cultivated land; in the area transferred from other land use types, the proportion of forest accounts for 49.59%. A total of 1825.45 km 2 of grassland was from other land use types, and 1198.30 km 2 was from cultivated land; the grassland transfer proportion is 35.57% of the whole basin. Although cultivated land was mostly

Calibration and validation of the SWAT model
The sensitivity analysis indicates that CN2 (SCS streamflow curve number for moisture Condition 2), CANMX (Maximum canopy storage), ALPHA_BNK (Baseflow alpha factor for bank storage), SOL_AWC (Available water capacity of the soil layer), SOL_K (Saturated hydraulic conductivity of the soil layer), SOL_BD (Moist bulk density of the soil layer), ESCO (Soil evaporation compensation factor), and REVAPMN (Threshold depth of water in the shallow aquifer for "revap" to occur) are more sensitive than the other parameters of the streamflow simulation. The sensitive parameters of streamflow in this paper are in accord with the summary for streamflow in the SWAT simulation of Athira (2021) [59], and their optimal values for the SWAT model are shown in Table 3. SWAT-CUP is used to calibrate and validate the model in this study; 1980-1985 was selected as the warm-up period, 1986-1994 was the calibration period, and 1995-1997 was the validation period. R 2 (the coefficient of determination), NSE (the Nash-Sutcliffe efficiency coefficients), BIAS (percent bias), P-factor and R-factor are used to assess the accuracy of the model simulation. When the R 2 and NSE values are more than 0.5, the BIAS value is less than or equal to ±20%, the P-factor is more than 0.6 and the R-factor is less than 1.5, the SWAT calibration results on a monthly scale are considered acceptable [60]. Table 4 shows the calculation results of the model evaluation. The goodness-of-fit statistics indicate reasonable agreement between the observed and simulated streamflows. However, there was a higher deviation in July in 1989 and 1996 because the precipitation on the 16th of July in 1989 was 26.7 mm and the precipitation on the 12th of July in 1996 was 91.9 mm, which accounted for 49.7% and 98.9% of the July precipitation of those two years, respectively, which may be the reason for the comparatively higher deviations between the observed and simulated streamflows. The SWAT model cannot accurately simulate rainstorm processes, and the heavy rains of 1989 and 1996 led to the poor simulation of those two years and reduced the overall accuracy of the simulation and validation periods.  , were less than those of other years, which is the same as the temporal characteristics of precipitation and indicating that the amounts of blue water were strongly related to rainfall, and the correlation coefficient of the two elements was 0.97. There is no significant relation between the amounts of rainfall and green water.
Blue water increases at rates of 0.87 million m 3 /yr and 2.2 million m 3 /yr in scenarios I and II, respectively; the green water increase rates are 10.29 million m 3 /yr and 10.8 million m 3 /yr in scenarios I and II, respectively. The increase rate of blue water in scenario II is approximately 2.52 times that of scenario I, but there is almost no difference in the increase rate of green water between the two different scenarios. Both blue water and green water show increasing trends in the two land use scenarios.
The blue water, green water and total water resources in scenario II were less than those in scenario I during years 38, 14 and 34 years of the 38 years studied (Fig 5c). The average annual differences in blue, green and total water resources between scenario II and scenario I are -72.08 million m 3 , 24.34 million m 3 , and -47.74 million m 3 , respectively, which indicates that land use change caused by the GGP leads to a decrease in blue and generalized water resources and an increase in green water resources.

• Interannual variations of the diverse hydrological components
To clarify the impact of land use on regional water resources, it was necessary to analyze and quantify the diverse components of the hydrological elements within the study area. They include SURQ, LATQ, GWQ, ET, and SW obtained from the well-calibrated SWAT model. Fig 6 demonstrates that the annual SURQ and component of blue water in scenario II were smaller than those in scenario I in all of the studied years except 2017; the annual average difference between scenario II and scenario I is 17.56 mm, and the maximum difference of 45.61 mm appeared in 2013. The change in SURQ was consistent with the fact that surface runoff has decreased on the Loess Plateau since the implementation of the GGP. LATQ in scenario II was larger than that of scenario I in 2000 and 2015 but lesser in the other years; the change range of LATQ was the smallest of the five variables, the annual average difference between scenario II and scenario I was 2.5 mm. The GWQ in scenario II was larger than that in scenario I in 1992,1995,2000,2005,2008 and 2015 and smaller in the other years; the annual average difference between scenario II and scenario I was 6.57 mm. The annual green water flow component of the green water in scenario II was larger than that of scenario I in all of the studied years except 2000, whereas the green water storage in scenario II was less than that of scenario I in all of the studied years (Fig 6). The mean annual green water flow had variations from 255.15 mm to 356.12 mm in scenario I and from 262.81 mm to 365.19 mm in scenario II; The mean annual green water storage had variations from 43.68 mm to 115.92 mm in scenario I and from 32.96 mm to 109.92 mm in scenario II throughout the study period; green water flows in scenario I were smaller than those of scenario II except in 1998, whereas green water storage values in scenario I were more than those of scenario II except in 1998. The differences of the green water flows between the two land use scenarios were comparatively larger than those of green water storage values. The annual average differences in green water flow and green water storage between scenario II and scenario I were 12.28 mm and -8.9 mm, respectively. The GGP reduced SURQ and green water storage but increases GWQ, LATQ and green water flow.
Spatial distribution changes of the blue and green water. Blue water, green water storage, and green water flow had considerable spatial variations among subbasins in the two land use scenarios (Fig 7). Water resource values were divided into five ranks according to the natural breaks method in ArcGIS software in Scenario I, whereas the classification boundaries in Scenario II were the same as those used in Scenario I to compare the spatial variations of the two scenarios.
The spatial distributions of the mean annual blue water were less than 59.32 mm, which is exactly the same in the two scenarios, and the ten minimum change regions were distributed in the upper area of the basin (Fig 7a). It can be observed that blue water in the second rank jumps from 59.33 mm to 209.72 mm, 14 subbasins in Scenario I and 36 subbasins in Scenario II were in that range. Thirty-seven subbasins have average annual blue water values ranging from 207.93 mm to 241.71 mm, while only 15 subbasins were in that interval in Scenario II. Fig 7a shows that blue water in Scenario II was less than that in Scenario I in all subbasins, which demonstrates that blue water decreased over the whole basin with the implementation of the GGP. The minimum changes occurred in the upper part, and the maximum changes occurred in the middle part of the basin. Fig 7b presents the green water storage across subbasins for the two scenarios. It indicates that the green water storage values lower than 47.76 mm were located in the upper reaches of the river basin in both of the two land use scenarios. The areas with values between 78 mm and 88.45 mm occupied 57.44 percent of the whole basin in Scenario II. Green water storage values in 8 subbasins in Scenario II were larger than those of Scenario I, and the area occupies 12.89 percent of the whole basin. The maximum decrease region is located in the northern and eastern parts of the basin, and the increase region is located in the western part of the basin.
The maximum green water flow values were located in the upper reaches in Scenario I but in the upper and middle reaches in Scenario II (Fig 7c). Data in 42 subbasins were more than 314.46 mm in Scenario II and the area accounted for 63.21 percent of the whole basin; green water flows lower than 276.64 mm accounted for 14.75 percent in Scenario II and 21.5 percent in Scenario I, indicating that there was comparatively little variability in areas of the minimum data range between two land use scenarios. Green water in one subbasin decreased, and its area was 0.79 percent of the basin; therefore, it can be regarded that the green water flow increased over the whole basin because of the GGP.
The land use conversion map from 1980 to 2017 and subwatersheds in the SWAT model are superimposed on one map (Fig 8) to compare land use conversion with spatial changes in water resources. We can see that spatial changes in blue water, green water storage, and green water flow (Fig 7) are related to land use conversion to a certain extent. The areas with the lowest reduction of blue water were consistent with the areas in which land use remained

PLOS ONE
Impact of the Grain for Green Project on water resources unchanged, such as the subbasins from one to ten, whereas the areas with the largest reduction in blue water were consistent with the areas in which land use changed significantly, such as subbbasins 17, 18, 19, 27, 28, 29, 32, 33, 35, 36, 38, where large amounts of cultivated land were converted into grassland and forestland. The second reduction areas, such as subbasins 11, 13, 16, 20, 21, and 23, were also the areas in which many cultivated lands were converted into grassland and forestland. The relationship between the spatial changes in green water flows and land use conversions were relatively the same as those between blue water values and land use conversions. There was no significant relationship between land use conversions and spatial changes in green water storage.

WF
Agricultural WF. This study calculated crop WF during the period of 1994-2017 because of crop output was inaccessible before 1994 (Fig 9). The agricultural WF in the Yanhe River Basin showed a rapid downward trend, with the highest value of 250 million m 3

PLOS ONE
Impact of the Grain for Green Project on water resources Table 5. The forest WF was larger in May, June and July because the actual evapotranspiration rates of these three months were greater than those of the other months, and the minimum value appeared in October because leaves begin to wither and fall, and the water demand is the smallest in the growth stage. The total forest WF increased from 5.22 billion m 3 to 19.58 billion m 3 during 1980-2017. After 2000, the forest WF began to increase rapidly, and the maximum value appeared in 2017.
The grass WF was larger in May and June but smallest in October in the grass growth stage, which is the same as the monthly characteristics of forests. Table 1 demonstrates that grassland area decreased from 3565.2 km 2 to 3547.09 km 2 , which indicates that GGP had comparatively little effect on grassland, so the grass WF had no obvious change and decreased slightly from 18.71 billion m 3 to 15.44 billion m 3 (Table 6).
Agricultural virtual water flow. The population in the Yanhe River basin increases from 60.92 to 79.83 ten thousand, but per capita grain consumption decreased from 208.04 kg to 108.5 kg, and the regional total grain consumption WF had a downward trend during 1994-2017. Crop yields varied greatly from year to year, leading to an irregular decrease in grain production WF. Fig 10 indicates that from 1994 to 2017, the grain consumption WF in the Yanhe River Basin ranged from 1.01 billion m 3 to 3.24 billion m 3 , and the grain production WF ranged from 1.36 billion m 3 to 2.5 billion m 3 . The grain production WF before 2000 was more than 2 billion m 3 and was approximately 1.5 billion m 3 after 20033 after . In 19943 after -19973 after , 20023 after -20033 after , 2012 and 2014-2017, virtual water flowed from the Yanhe River basin to the outside. From 1994 to 2017, approximately 8 million m 3 of virtual water came from outside the region annually. Therefore, virtual grain water in the Yanhe River Basin imported from the outside, which can alleviate the pressure of water resources in the basin to a certain extent.  Table 7. The total WF showed an increasing trend due to the significant increase in forest WF, although agricultural and grass WFs showed decreasing trends during the study years. However, the ecological water stress index had no obvious temporal change because the change trend of generalized water resources was the same as that of the regional total WF. The EWSI was less than 1 in the five studied years in both land use scenarios from the point view of generalized water resources. The ESWI in scenario II was larger than that in scenario I in the five  years, and this demonstrates that the GGP slightly increased regional water stress because the difference in ESWI is approximately 0.01 between the two land use scenarios. Table 4 and Fig 4 indicated a good SWAT model performance for streamflow simulation in the Yanhe River basin. However, factors such as parameter selection, rainstorm, observation data missing and heterogeneity within HRUs can all lead to deviation of simulated data from observed data. Currently, SWAT Calibration and Uncertainty Programs (SWAT-CUP) is mostly to determine the sensitive parameters and the optimal value. CN2 is mostly recognized as the most important parameter in SWAT streamflow simulations [61,62] because SWAT uses the soil conservation service streamflow curve method to calculate the surface streamflow. In addition, the surface streamflow in most watersheds plays a dominant role in the total annual streamflow. Therefore, CN2 plays important role in streamflow simulation and it is critical for model accuracy. Previous research showed that there was a great difference in parameter of CN2 in different study in Yanhe river basin (e.g. Studies have shown that multiple outlets in SWAT model simulation perform better than a single outlet [75,76]. Ganguyi hydrology station was a single calibration outlet in this paper that may be one of the reasons for relative low accuracy. Simulation accuracy of areas for severe land use changes is smaller than that of areas with gentle land use change because land use cannot be input into the model every year. For example, simulation results in Xihe River basin [77] and Qingshuihe basin [78] were obviously better than that in the Loess Plateau, The upper reaches of the simulation results was better than the middle and lower reaches in the Yellow River basin [79].

Assessment of the water balance components and uncertainty analysis
This study investigated the spatial-temporal impact of the GGP on regional water resources and ecological water stress in the Yanhe River Basin in 1980 and 2017 using the results of the SWAT model. Many researchers have explored blue water and green water resources since the appearance of the concept using different hydrological cycle models [43]. Rost et al. [80] used the Lund-Potsdam-Jena managed land (LPJmL) model to quantify the global consumption of both blue and green water in 1971-2000. SWAT is currently the most commonly used model to simulate the effects of land-use change and climate change on the hydrological cycle. Many previous studies have demonstrated that the results of SWAT simulations are reliable when exploring land-use and climate change effects on the hydrological cycle in many river basins around the world [81,82].
The most important elements of the water balance for one watershed include precipitation, SURQ, LATQ, GWQ, ET and SW [83]. In this study, hydrological components other than precipitation were analyzed to investigate hydrological cycle changes caused by the GGP. Average annual water balance components simulated by the SWAT model is presented in Table 8. This result shows that SRUQ decreased, ET increases were accompanied by vegetation cover and urban land use increases in the Yanhe River basin, which is consistent with the same phenomenon found in many regions around the world [84][85][86]. ET was the largest contributor to water loss in this region. There was only an 8.49 mm difference in water yield between Scenario II and Scenario I because SURQ decreased while GWQ and LATQ increased after the GGP. The drastic interannual changes of water yield in both scenarios may be caused by parameters other than land use change.
In this study, the accuracy of SWAT simulation was one of the important key factors because many results, such as blue water, green water, and hydrological elements (green water flow, green water storage, SURQ, GWQ and LATQ), were all outputs of the model. The R 2 and NSE values in this study were within the accuracy requirement of model establishment, and the accuracy was close to the results of other studies in this river basin [87][88][89]. During the SWAT simulation, land use was classified into six types that cannot actually represent the land use type, and the land use module should be improved. Data from four meteorological stations, Yanchang, Ansai, Zhidan and Yanan, were input into SWAT, and the precipitation data was also from meteorological stations. More meteorological and precipitation data may enhance the simulation accuracy. Reservoir information, such as Wangyao, can be used in the model to consider human impact, while improving DEM data can also improve simulation accuracy [90].
Only five crops, wheat, summer maize, potato, soybean and oil crops, were considered, while other agricultural WFs were not included because of the lower yield in the Yanhe River Basin; therefore, the agricultural WF was less than the actual total agricultural WF in the region. The most important varieties of trees in the Loess Plateau are robin pseudoacacia, Platycladus orientalis, Pinus tabulaeformis and so on. Different experiments on forests on the Loess Plateau show that the ecological water demand coefficient of arbor forests is 0.757 and that of shrub forests is 0.612 [91]; the water demand coefficient of forestland is 0.76, that of shrub forests is 0.61 and that of sparse forestland is 0.48 [92,93]. There are no precise and uniform water demand coefficients of forest and grass on the Loess Plateau because there have been few experimental studies on it. This paper cites the existing study results and does not distinguish forest types because of a lack of experiments and image interpretation accuracy.
The above factors will all influence the ecological WF and EWSI.

Implication of the GGP on water resources system
The results of this study show that blue water resources decreased and green water increased obviously, but generalized water resources decreased slightly; SURQ and SW decreased but ET increased significantly over the whole basin with the achievements of the GGP; these results are similar to those of Yang [87] and Yang [94] in the Yanhe River Basin. Large amounts of cultivated lands were transformed into forestlands, reducing slope surface runoff and ultimately improved green water storage, which is important for restoring vegetation within this region [95]. The reason that the blue water decreased and green water increased may be the shift from the former to the latter. With the increase in vegetation coverage in the Yanhe River Basin, the accumulation of bark debris and leaf litter significantly increased soil surface roughness, slowing the runoff rate, increasing infiltration and intensifying evaporation [96]. The water demand of increasing vegetation in the semiarid climate of the Loess Plateau has reduced the soil water content of the top 3-5 m to almost withering humidity and further exacerbated soil drying [97]. The ET in the growing season limits vegetation growth [98]. Water resources cannot satisfy water requirements due to forest and grass growth leading to slow growth of trees and low yields of forests because precipitation in this region did not significantly increase during the study period although evaporation increased. The Chinese government plans to invest another US$9.5 billion in the GGP on the Loess Plateau by 2050 [99]. The results of this study show that although this project has contributed to increasing forestland https://doi.org/10.1371/journal.pone.0259611.g011 area and green water resources, this has been at the cost of detectable reductions in blue water resources, especially in river runoff. Studies have indicated that current vegetation productivity in the Loess Plateau is already close to net primary productivity (NPP), and total water resources will inevitably decrease for human use to less than the required amount if vegetation coverage continues to increase [22]. Therefore, future GGPs should consider water stress and other negative impacts.

GGP sustainability
Natural rainfall data from 43 articles and 331 runoff experimental graphics in the Loess Plateau have been integrated and synthesized to analyze the land use change impact on runoff, and the results show that artificially promoted vegetation has a substantial negative impact on runoff formation [100]. Blue water decreased in all of the studied years, and the difference between the two land use scenarios is -72.08 million m 3 . Vegetation cover increases from 4347.6 km 2 to 6565.61 km 2 . SURQ decline reduced the amount of water resources that people can use directly; hence, during the implementation of the GGP more detailed investigation about how to adjust vegetation types according to regional precipitation, topographical conditions and temperature should be carried out. The average green water increase indicates that the ecological environment has improved significantly, and the EWSI became slightly worse after the GGP. Therefore, vegetation restoration has both positive and negative impacts on water resources and ecological systems. Appropriate trees or other vegetation types should be planted to maintain the balance between ecosystems and water resource systems in the future.
Compared with the previous research on impact of land use change on regional water resources, this paper evaluates the role of green water resources in the basin; analyze regional water stress variation adopting the concept of generalized water resources. Taking green water as one part of regional water resources and taking forest and grass land water footprint as one part of water requirement, this method can discuss whether vegetation restoration aggravate or alleviate regional water resource stress more comparatively.

Conclusions
The concern about the impact of vegetation restoration on regional water resources has increased globally in recent years, and it is necessary to integrate the SWAT model and WF to explore the problem. This paper presents a framework to analyze the spatial and temporal characteristics of blue water, green water and each hydrological element as well as EWSI changes. The Yanhe River basin is taken as a case study to explore the negative and positive impacts of the GGP on regional water resources. The results show that the vegetation cover area increased significantly since the implementation of the GGP, blue water, and SURQ decreased substantially, which is consistent with the results of many previous investigations. Green water and generalized water resources had increasing and decreasing trends, respectively, which were caused by land use change during 1980-2017. The EWSI became slightly worse from the perspective of the generalized water resources and the WF concept.